# The following script was adopted from Thiele, Jan C., Kurth, Winfried and Grimm, Volker (2014) 
# 'Facilitating Parameter Estimation and Sensitivity Analysis of Agent-Based Models: 
# A Cookbook Using NetLogo and 'R'' Journal of Artificial Societies and Social Simulation 17 (3) 11 
# <http://jasss.soc.surrey.ac.uk/17/3/11.html>.

# a function to handle a simulation

simulate <- function(param.set, parameter.names, no.repeated.sim, nl.obj, trace.progress, iter.length, function.name) {
  # some security checks
  if (length(param.set) != length(parameter.names))
  { stop("Wrong length of param.set!") }
  if (no.repeated.sim <= 0)
  { stop("Number of repetitions must be > 0!") }
  if (length(parameter.names) <= 0)
  { stop("Length of parameter.names must be > 0!") }
  
  # an empty list to save the simulation results
  eval.values <- NULL
  
  # repeated simulations (to control stochasticity)
  for (i in 1:no.repeated.sim)
  {
    # create a random-seed for NetLogo from R, based on min/max of NetLogo's random seed
    
    NLCommand("random-seed",runif(1,-2147483648,2147483647), nl.obj=nl.obj)
    
    NLCommand("setup", nl.obj=nl.obj)
    
    # set NetLogo parameters to current parameter values
    lapply(seq(1:length(parameter.names)), function(x) {NLCommand("set ",parameter.names[x], param.set[x], nl.obj=nl.obj)})
    
    # 24 months for start up
    
    NLDoCommand(1,"repeat 12 [go]", nl.obj=nl.obj)
    
    # 144 months for 12 years model simulation (12 x 12)
    
    cal.crit <- NLDoReport(12,"repeat 12 [go]",c("unemployment-rate","ltu-rate", "theta", "wag", "ltu-wag", "dur", "ltu-dur"), 
                           as.data.frame=T, df.col.names=c("urate","lturate", "theta", "w", "ltuw", "d", "ltud"), nl.obj=nl.obj)
    
    # calculate calibration criteria
    
    umean.criterion <- mean(cal.crit$u) 
    theta.criterion <- ( sd(cal.crit$theta) / mean(cal.crit$theta) )
    ltu.criterion <- mean(cal.crit$lturate)
    wage.criterion <- mean(cal.crit$w)
    ltu.wage.criterion <- mean(cal.crit$ltuw)
    duration.criterion <- mean(cal.crit$d)
    ltu.duration.criterion <- mean(cal.crit$ltud)
    
    # merge calibration criteria
    calibration.criteria <- c(umean.criterion, theta.criterion, ltu.criterion, wage.criterion, ltu.wage.criterion, 
                              duration.criterion, ltu.duration.criterion)
    
    # append to former results
    eval.values <- rbind(eval.values,calibration.criteria)
  }
  
  # print the progress if requested
  if (trace.progress == TRUE)
  {
    already.processed <- get("already.processed",env=globalenv()) + 1
    assign("already.processed", already.processed, env=globalenv())
    print(paste("processed (",function.name,"): ", already.processed / iter.length * 100, "%", sep = ""))
  }
  
  # return the mean of the repeated simulation results
  if (no.repeated.sim > 1)
  {
    return(colMeans(eval.values))
  }
  else {
    return(eval.values)
  }
}
